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In transport calculations for molecular junctions based on density functional theory the choice of exchange 
and correlation functional may dramatically affect the results. In particular local and semi-local functional 
tend to over-delocalize the molecular levels thus artificially increasing their broadening. In addition the same 
molecular levels are usually misplaced with respect to the Fermi level of the electrodes. These shortfalls are 
reminiscent of the inability of local functionals to describe Mott-Hubbard insulators, but they can be corrected 
with a simple and computationally undemanding self-interaction correction scheme. We apply such a scheme, 
as implemented in our transport code Smeagol, to a variety of phenyl-based molecular junctions attached to gold 
electrodes. In general the corrections reduce the current, since the resonant Kohn-Sham states of the molecule 
are shifted away from the contact Fermi level. In contrast, when the junction is already described as insulating 
by local exchange and correlation potentials, the corrections are minimal and the I-V is only weakly modified. 



PACS numbers: 



o 



> 

O 



X 



I. INTRODUCTION 

Devices in which organic molecules form the active ele- 
ment are taking an increasingly important place in solid state 
physics. These are the building blocks of the field of molec- 
ular electronics, whose potential applications underpin a vast 
number of technological areas. Novel components for high- 
performance computers', highly sensitive chemical sensors^, 
and disposable and wearable devices are among the many pro- 
posed applications. In addition the same set-up can be used for 
medical purposes, for instance in the detection of viruses^ or 
in the construction of a rapid and reUable protocol for DNA 
sequencing^. 

Most of the progress originates from the recently achieved 
ability to construct single molecule junctions and to measure 
their transport properties. The experimental strategies for fab- 
ricating such devices include mechanically controllable break 
junctions^'*"'^, scanning tunneling spectroscopy^ *^, lithograph- 
ically fabricated nanoelectrodes^ and colloid solutions'". Un- 
fortunately most of these techniques are actually "blind" in 
the sense that the geometrical configuration of the device at 
the atomic scale is not known. For this reason it should 
not be surprising that different experiments for the same 
molecule yield conductances differing by of up to three or- 
ders of magnitude^'^'** '". In view of all these uncertainties ab 
initio quantum transport schemes for calculating the electrical 
response of a molecular device become crucial. 

The most common computational scheme for evaluating 
the electronic transport of a molecular device combines scat- 
tering theory in the form of the non-equilibrium Green's 
function (NEGF) formalismii, with an electronic structure 
method, the most widely used being density-functional the- 
ory (DFT)'^ '^ '"*. Alternatives are time-dependent DFT'^ or 
many-body methods'^, which at present however are more 
computationally involved and are at an earlier stage of de- 
velopment as routinely used schemes. Unfortunately for sev- 
eral molecules NEGF combined with DFT appears to disagree 
with the experimental results. In particular it seems to system- 



atically overestimate the current flowing across a device, even 
if one bears in mind the uncertainty over the detailed geome- 
try of the junction. In the prototypical case of benzenedithiol 
(BDT) sandwiched between two gold electrodes, the conduc- 
tance for the most probable contact geometry calculated with 
NEGF and DFT is higher than that of any of the experiments 
by at least one order of magnitude '^•'**'"'^-. 

In this work we demonstrate that most of the eiTors origi- 
nate from using local or semi-local exchange and correlation 
potentials in DFT. This shortfall however can be coiTected 
by applying an atomic self-interaction correction (ASIC)^ 
scheme to the local functionals, leading to a dramatic im- 
provement of the agreement between the calculated results 
and those obtained from experiments. The reason for this 
improvement is rooted in the ability of ASIC to predict the 
coiTect ionization potential of molecules, and hence to cor- 
rectly describe the band alignment between the molecular or- 
bitals and the Fermi level, Ep, of the metallic electrodes. 
Here we extend the results previously published for BDT 
attached to the gold hollow site of the (111) surface^^ to 
new molecules such as benzenedimethanethiol (BDMT) and 
biphenyldithiol (BPD). We also carry out a thorough inves- 
tigation of the effects of changing the anchoring geometry 
between the molecule and the electrodes, demonstrating that, 
for the case of BDT attached to gold, the ASIC calculated I- 
V curves are rather stable with respect to geometry changes. 
This explains the relatively narrow distribution of conduc- 
tance histograms found in experiments^. 



II. METHOD 

A. Non-Equilibrium Green's Function Formalism 

The problem we wish to solve is that of calculating the two- 
probe I-V curve of a molecule sandwiched in between two 
metallic electrodes. The NEGF scheme partitions such sys- 
tem into three regions, respectively the two current/voltage 
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electrodes (leads) and a middle region called the scattering re- 
gion (SR) (see figure[T]i. 
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FIG. 1: Schematic diagram of a metal-molecule-metal junction. A 
scattering region is sandwiched between two current/voltage probes 
kept at the chemical potentials /ii and ^2- The electrodes are mod- 
eled as being periodic in the direction of transport. A number of 
layers of the electrodes are included in the scattering region (the part 
of the system enclosed by the dashed box) to allow the charge density 
to converge to the bulk value. 

The SR includes the molecule and a portion of the leads, 
which is extended enough for the charge density calculated 
at the most external atomic layer to resemble that of the bulk 
electrodes. The leads, assumed to be a periodic crystal in the 
direction of the transport, are kept at different chemical poten- 
tials /ii/2 = E-p ± eV/2 where V is the applied potential bias 
and e is the electron charge. The SR is described by a Hamil- 
tonian Hs, which is used to construct the non-equilibrium 
Green's function, G{E) 

G{E) = \im[{E + iT^)-H,-^i-Y.2]-\ (1) 

where S1/2 are the self-energies for the leads. G{E) enters 
in a self-consistent procedure to calculate the density ma- 
trix, p, of the SR, and hence the two probe-current of the 
devic o"''^''^-^^ . The non-equilibrium p is calculated follow- 
ing the NEGF prescription as 

1 1"°° 

P^^J G{E)[Tif{E,fii)+T2f{E,p2)]G\E)dE, 

(2) 

where ri/2 ~ i[^i/2 — ^1/2]- practice, this integral is per- 
formed by splitting it into two part s"''-^''"^'^^ : an equilibrium 
part which can be integrated along a contour in the complex 
plane, and a non-equilibrium part which has to be integrated 
along the real energy axis but which only contributes around 
Ep. The non-equilibrium charge density p is then used to 
calculate a new Hamiltonian for the scattering region Hs[p] 
(where it is assumed that the Hamiltonian has some func- 
tional dependence on p). This procedure is repeated self- 
consistently until the density matrix converges. Finally, the 
converged Green's function can be used to calculate the cur- 
rent / through the device: 

Op r°° 

/ = - y Tr[G{E)TiGHE)T2]{fiE, fii)-f{E, ii2))dE. 

(3) 



This is effectively the integral between p,i and p,2 (the 
bias window) of the transmission coefficients T{E) — 
Tr[G{E)TlG^E)T2]. In brief the NEGF scheme calculates 
the non-equilibrium scattering potential of the device. Thus 
the transmission coefficient T{E) is simply a superposition 
of resonances located at the molecular single-particle energy 
levels, including a possible shift and broadening due to the in- 
teraction with the leads. It follows that the molecular levels 
close to the electrode Fermi level, i.e. the highest highest oc- 
cupied molecular orbital (HOMO) and the lowest unoccupied 
molecular orbital (LUMO), provide the dominant contribution 
to the current. It is therefore crucial to employ an electronic 
structure theory capable of describing the position of these 
two levels accurately. 
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FIG. 2: Energy level diagram of the metal-molecule-metal junction 
shown in figure [T] (a) Energy level line up and (b) transmission co- 
efficients as a function of energy. The resonances in the transmission 
coefficients correspond to energy levels in the molecule. 



B. DFT and the self-interaction problem 

The NEGF scheme is general and not related to a particular 
electronic structure theory. However, since the transmission 
occurs through single-particle states, the associated electronic 
structure theory should meet several requirements^''. First of 
all, the single-particle levels must closely resemble the phys- 
ical removal energies of the system, i.e. the energy levels 
of the molecule should line up correctly with the Ep of the 
metals forming the electrodes. Secondly, the theory should 
work at both integer and fractional occupation, ensuring the 
correct response of the molecular orbitals to changes in occu- 
pation due to the applied potential bias. Finally, the matrix 
elements describing the interaction between the leads and the 
molecule should be calculated accurately avoiding the use of 
phenomenological parameters or fitting procedures. 

The Kohn-Sham (KS) form^^ of DFT is certainly the most 
widely used ab initio electronic structure theory associated 
with the NEGF method. In the KS scheme the DFT prob- 
lem of finding the ground state charge density is mapped 
onto that of solving a system of non-interacting single-particle 
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Schrodinger equations 

[„iv2 + v{r) + u{p) + v^c{p, r)]C(r) = e^^,^(r) , (4) 

where u{p) is the Hartree (classical Coulomb) potential and 
vy^c is the exchange-correlation (XC) potential. The exact 
form of vxc is unknown, so it is typically approximated by 
local functionals such as the local density approximation^'' 
(LDA) or the generalized gradient approximation^'' (GGA). 

Therefore the KS scheme is only a convenient prescription 
for minimizing the energy functional, and as such the individ- 
ual KS eigenvalues do not necessarily correspond to physical 
energy levels of the system. The DFT observables such as 
the total energy and the charge density are in fact only inte- 
gral quantities respectively of the KS eigenvalues and eigen- 
vectors. However, in solid state theory it is common to asso- 
ciate the KS energies with the system band structure. This, 
although it is not justified at the fundamental level, is sup- 
ported by the practical evidence that KS bands represent a 
good first approximation of the true energies of a system, 
particularly for metals. However, there is a remarkable ex- 
ception to the non-physical nature of the KS eigenvalues. In 
fact, the KS energy of the HOMO level (e^QMo) tie rig- 
orously associated with the negative of the ionization poten- 
tial 7^^ -**. Similarly the HOMO KS energy for the negatively 
charged system can be associated with the chemical affinity A. 
This suggests that, at least for moderate bias where the trans- 
port is through the HOMO, KS theory can be used effectively 
for transport calculations^^. Unfortunately, standard local and 
semi-local functionals completely misplace e^oMO' which is 
often nowhere near — /. This is particularly problematic for 
organic molecules, for which the LDA Chomo i^ typically off 
from the experimental value of / by approximately 4 eV^'. 

Most of the failures of LDA and GGA can be traced back to 
the self-interaction (SI) problem, i.e. the spurious interaction 
of an electron with the Hartree and exchange-correlation (XC) 
potentials generated by itself'"'. In the case of Hartree-Fock 
theory the self-Hartree energy U[p'^ cancels exactly the self- 
XC energy ^xc[pj^,0] 

i?xcK,0] + C/K]-0, (5) 

where we introduce the orbital density pj^ — I'f/'J^p. How- 
ever, this cancellation is incomplete for both LDA and GGA. 
The resulting KS potential therefore appears too repulsive and 
the eigenstates of a molecule are pushed to higher energies. 
When this is transferred to the transport problem, the posi- 
tion of the peaks in the transmission coefficient results are in 
the wrong place. In particular it is likely that peaks will erro- 
neously move close to the Fermi level of the electrodes, gen- 
erating a large conductivity at zero bias. Since the position of 
the peaks is ultimately determined by the interplay between 
the hopping probabilities from the electrodes to the molecule 
with the charging energy of the molecule itself, this failure is 
somewhat reminiscent of the inability of local and semi-local 
functionals to describe Mott-Hubbard insulators. 

Perdew and Zunger suggested^" that the self-interaction 
correction (SIC) 5^^'^ for an occupied KS orbital ij}'^ can be 



calculated from the sum of the self-Hartree and self-XC ener- 
gies 

5T^E]i^^[plM + U[pl]. (6) 

Note that, although here we consider the case of SIC-LDA 
only, the same procedure can be readily applied to any ap- 
proximated functionals. The SIC-LDA XC energy is 
thus obtained by subtraction 

occ. 

EfS[p'.p']=Ek^\p\p']-Y.^f''^ il) 

n<7 

and the SIC-LDA XC potential w^^'^ '^, to be subtracted from 

LDA.CT .... , 

t^xc 7 IS Simply given by 

«f/c- = «([p„];r)+4g^-'^([pT,0];r), (8) 

with 

ui[p];r) = /dV^, (9) 

z;^°^-([pT,pi];r) = ^^Ek^^[p\pi] . (10) 

The problem of finding the energy minimum is complicated 
by the fact that E^^ is not invariant under unitary rotations of 
the occupied KS orbitals, which instead leave p invariant. This 
creates problem to the standard KS scheme since the theory 
becomes size-inconsistent. In order to avoid such a complica- 
tion modern SIC theory introduces a second set of orbitals (p^^ 
related to the canonical KS orbitals tp^ by a unitary transfor- 
mation M 

m 

The functional can then be minimized by varying both the KS 
orbitals and the unitary transformation M, leading to the sys- 
tem of KS-like equations: 

= [HS + Avf^'mr) = 6;^-Sic^-(r) , (12) 

where Hq is the LDA Hamiltonian, and the SIC potential, 
AwfJ^, is given by: 

A.^ic = Y: M^^^vf^^^ = E vf^^P:t , (13) 

where is the projector |05^)(0f„|. 

The numerical implementation of the full Perdew-Zunger 
SIC scheme in typical solid state codes is cumbersome since 
the theory is orbital dependent and the energy minimization 
cannot follow the standard KS scheme^'. When this is ap- 
plied to transport there is the additional complication that the 
KS orbitals are never individually available and moreover that 
one has to deal with an intrinsically non-local potential. A 
drastic simplification of the problem can be obtained by us- 
ing a recently implemented atomic approximation to the SIC 
scheme, which we call the ASIC^'. This is based on the 
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pseudo-SIC approximation, originally proposed by Vogel and 
coworkers^, and later extended to non-integer occupation by 
Filippetti and Spaldin^. The main idea is that of replac- 
ing the (pm orbitals with atomic like functions, which are not 
calculated self-consistently. Thus the SI correction becomes 
atomic-like and no information is needed other than the charge 
density and the ASIC projectors^. 

ASIC has been demonstrated to produce single-particle 
energy levels which match the experimental molecular re- 
moval energies quite wel^. In particular, for several differ- 
ent molecules investigated, ej^Qj^Q w — /, and the Ehomo 
for negatively charged molecules is close to the molecular 
affinity^'. As an example, ASIC places Ehomo ^'^^ 1,2-BDT 
at -8.47 eV to compare with the LDA value of -4.89 eV and 
the experimental ionization potential ^8.5 eV^. ASIC is the 
scheme that we adopt in this work. 



C. ASIC and the energy derivative discontinuity 



D. Calculation Details 

We have numerically implemented the ASIC method^l in 
the localized atomic orbital code SIESTA^'*, which is the DFT 
platform for our transport code Smeagol^^, and carried out cal- 
culations for the prototypical Au/Benzene/Au molecular de- 
vices. The ASIC corrections are not applied to the gold atoms, 
as SI errors for metals are small^'. Unless otherwise speci- 
fied, we use a double zeta polarized basis set^"* for carbon and 
sulfur s and p orbitals, double zeta for the Is orbitals of hy- 
drogen and 6s-only double zeta for gold. The mesh cut-off 
is 200 Ry and we consider 500 real and 80 complex energy 
points for integrating the Green's function. Calculations were 
also performed using double zeta 6 s and single zeta 5d and 6p 
orbitals for the gold to investigate the effect of an enlarged ba- 
sis on the transport. Results for calculations using both basis 
sets are shown, with the 6s-only basis used unless otherwise 
indicated. For geometry optimizations and total energy cal- 
culations, the extended 5d6s6p basis set is employed, as the 
6s-only does not yield enough accuracy. In calculating the I- 
V curves, the potential bias was incremented in steps of 0. 1 
Volts. 



Local and semi-local functionals are affected by another 
fundamental problem, i.e. the lack of the derivative discon- 
tinuity (DD) in the DFT energy. This is the discontinuity at 
integer occupation in the derivative of the total energy of a sys- 
tem E{N) with respect to its occupation N—. From Janak's 
theorem it follows that there should be also a discontinuity in 
^HOMO (^) when going from N to (N + S) with 6^0+. 

In a previous work^"*, we have investigated the conse- 
quences of the absence of the DD in local XC functionals 
on the electronic transport of organic molecules. We have 
demonstrated that the DD can strongly affect the I-V char- 
acteristics of metal-molecule-metal junctions when the cou- 
pling between the molecule and the metal is relatively weak. 
This would be the case, for instance, where molecules bind to 
adatoms on the metal surface. However, we have also shown 
that in the case of strong coupling the DD had httle effect on 
the transport. 

The self-interaction error is largely responsible for the dis- 
appearance of the DD in local and semi-local functionals. It 
follows that SIC methods restore the DD at least in part. This 
last aspect is deeply rooted in the orbital minimization needed 
to extract the SIC potential, which unfortunately is not cap- 
tured by our ASIC approximation. Therefore the main feature 
of our ASIC methodology is simply the correction of the ion- 
ization potentials of the molecules^', thus yielding a quantita- 
tively realistic energy level alignment. Hence, the calculations 
described in this paper investigate a second important aspect 
of the SIC beside that of the DD. This is the quantitative de- 
scription of the band alignment between the metal Fermi en- 
ergy and the molecular levels. In particular, for most of the 
paper we investigate molecules attached to the gold fee (111) 
hollow sites, for which the coupling is expected to be strong 
and the effects due to the DD small. 



III. RESULTS 

We compare the electronic transport properties for 
three different molecules calculated either by using LDA 
or ASIC. These ai-e 1 ,4-benzenedithiol (1,4-BDT), ben- 
zenedimethanethiol (BDMT), and biphenyldithiol (BPD). 
The electronic transport through these molecules have al- 
ready been investigated at length both experimentally and 
computationally^i^'i^i^i'^' '°i '^1^'^- ''^1^°. Thus these offer us a 
unique benchmark for our calculations. In all cases, the 
molecule is attached via the sulfur atoms to fee (111) gold 
electrodes on each side. In the case of BDT, we look at sev- 
eral different anchoring geometries. It is worth pointing out 
that although the exact anchoring geometry encountered in 
breaking junction experiments is unknown, the resulting con- 
ductance histograms are relatively narrow suggesting an in- 
trinsic stability in the measurements^. This may be an indica- 
tion that similar anchoring configuration are highly probable 
in the breaking process, or alternatively that several anchor- 
ing geometries yield the same I-V characteristic. One of our 
goals is to distinguish between these two possibilities. 



A. Benzenedithiol 

The first device which we consider is 1 ,4-benzenedithiol 
(1,4-BDT) attached to two fee (111) gold surfaces. The an- 
choring geometries for the sulfur atoms investigated are the 
hollow site, the top site, the bridge site, (see figure O, as well 
as asymmetric configurations where the S ions are attached 
to an adatom on one side and to the hollow site on the other. 
We also examine the effect of altering the angle which the 
molecule makes with the metal surface, and of varying the 
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distance between the sulfur atom and the surface (i.e. vary- 
ing the strength of the coupHng between the molecule and the 
metal). Finally, we investigate the effect of hydrogenating the 
thiol groups. 



a) 
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FTG. 3: (Color on line). The different anchoring geometries investi- 
gated. Panel (a) defines the "sulfur-surface separation", dss, as the 
distance between the S anchoring ion and the Au fee (111) plane. 
In (b) we show the different possible anchoring sites. A is the hol- 
low site, B is the top site, and C is the bridge site. Color code: Au 
atoms=yellow (or light gray), S atoms=brown (or dark gray). 

As explained previously the actual experimental contact ge- 
ometry is unknown. Electronic structure calculations indicate 
that the lowest energy configuration for the molecule on the 
(111) surface occurs when the S atom attaches to the hol- 
low site, although other calculations suggest that the bridge 
site configuration has a lower energy Recent X-ray stand- 
ing wave measurements suggest that molecules in monolayers 
prefer to attach to adatoms on the metal surface'"'. It is also 
worth noticing that in breaking junctions the actual geometry 
might be different from any equilibrium geometries, since the 
system is likely to be under strain. Hence, a thorough analysis 
of this system should comprise many different configurations. 

We start by attaching BDT at the Au (111) hollow site 
(see figures [3] and nil. The distance of the sulfur atom from 
the plane of the gold surface (the "sulfur-surface separation", 
dss, shown in figure [3^) is optimized (using the 5d6s6p ba- 
sis set for gold) to a value of 1.9 A. This corresponds to a 
distance of 2.53 A between the sulfur atoms and the nearest 
gold atoms on the surface in good agreement with previous 
calculation s ' ^i^^i^^ . 




FIG. 4: (Color on line). BDT molecule attached to the hollow site 
of the Au (111) surface. The sulfur-surface distance, dss, is 1.9 A. 
Color code: Au=yellow (or light gray), C=black, S=dark yellow (or 
gray), H=blue (or dark gray). 

The I-V characteristic, the orbital resolved density of states 
(DOS) and the conductance as a function of energy, G, are 
presented in figure |5] The conductance is plotted in units of 



=j^, so that it is equivalent to the transmission coefficients 
T. LDA results are shown in the left panels, while the ASIC 
ones are on the right. From the DOS it is clear that the effect 
of ASIC is that of shifting the occupied orbitals downwards 
in energy relatively to the Ep of gold. The ASIC HOMO- 
LUMO gap is considerably larger than that calculated with 
LDA and most importantly in the case of ASIC there is lit- 
tle DOS originating from the molecule at E-p. This has pro- 
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FIG. 5: Transport properties of a BDT molecule attached to the gold 
fee (111) hollow site. The left plots correspond to LDA and the right 
ones to ASIC. The upper panels are the DOS of the S and C tt or- 
bitals ((a) and (b)), the middle are the transmission coefficients as a 
function of energy for various bias ((c) and (d)) and the lower are the 
I-V curves. Figure (f) is a zoom of (e) and compares our results with 
experiments from reference-. The vertical lines in (c) and (d) mark 
the bias window. 

found effects on the electron transmission. The LDA peaks of 
T{E) arising from occupied orbitals are shifted downwards 
in energy and away from E-p. At variance from LDA (figure 
|5};), where T{Ep) is dominated by a resonance at ehomo, the 
ASIC transmission (figure |5}l) is through the BDT gap and 
therefore it is tunneling-like. This results in a drastic reduc- 
tion of the low-bias current when going from LDA to ASIC 
(figure|5^). The ASIC-calculated conductance at zero bias is 
now about O.O6G0 (Go — 2e^//i), compared to 0.23Go for 
LDA. A conductance of O.O6G0 is much closer to the value of 
0.011 Go obtained by Xiao et. al.^ and is actually lower than 
values 0.09-0. 14Go obtained by Tsutsui et. ali^. 

Without altering the anchoring geometry, the basis set on 
the gold atoms is then changed to include the 5d and 6p or- 
bitals. The relative orbital resolved DOS, transmission coeffi- 
cients and I-V curves are presented in figure|6]for both LDA 
and ASIC. Clearly the enriched basis set does not have a large 
effect on the electronic transport, particularly at low bias. As 
can be seen from figure |6^, the I-V curves calculated with 
the 6s-only basis set for gold are approximately the same as 
that calculated with the 5d6s6p basis set up to about 1 Volt for 
both the LDA and ASIC. Differences appear only at high bias 
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(T^ > 1 Volt) and are due to the presence of Au d electrons 
in an energy range comprising the bias window. However the 
zero-bias conductances, from the transmission plots in panels 
(c) and (d) of figure |6l are very similar to those for the 6s- 
only basis set with a value of 0.21 Go for LDA and a value of 
O.O6G0 for ASIC. This demonstrates that the 6s-only basis set 
should give reasonably reliable results for electronic transport 
at low bias. 
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FIG. 6: Transport properties of a BDT molecule attached to the Au 
fee (111) hollow site calculated with an enriched 5d6s6p basis set for 
gold. The left plots correspond to LDA and the right ones to ASIC 
results. The upper panels are the DOS of the S and C n orbitals ((a) 
and (b)), the middle are the transmission coefficients as a function of 
energy for various bias ((c) and (d)) and the lower ((e) and (f)) are the 
I-V curves, including a comparison with the results shown in figure 
|4]for the 6s-only basis. Figure (f) is a zoom of (e) and compares our 
results with experiments from reference^. The vertical lines in (c) 
and (d) mark the bias window. 

Next, the transport properties of the system are calculated 
for the hollow site at different dss, as well as for the equi- 
librium distance c?ss = 1-9 A but different angles of the 
molecule with respect to the direction of transport. The re- 
sults are reported in figure |7] In general ASIC seems to be 
much less sensitive to changes in the anchoring geometry than 
LDA, particularly for the case of bond stretching. For in- 
stance, the zero bias conductance values calculated with LDA 
are 0.16 Go, 0.23 Go, 0.32 Go and 0.77 Go for dgs values of 
1.8 A, 1.9 A, 2.1 A and 2.5 A respectively. For the same dss 
ASIC returns 0.05 Go, 0.06 Go, 0.07 Go, and 0.14 Go- 

The stability of the ASIC I-V curve with respect the details 
of the anchoring geometry of the hollow site is an interesting 
feature, since it is suggesting that in breaking junction exper- 
iments several different anchoring configurations may yield 
similar I-V curves. This is consistent with the relatively nar- 
row peaks measured in the typical conductance histograms^. 
LDA does not seem to have this property and in fact the zero- 
bias conductance increases quite drastically when dss is in- 
creased. While counterintuitive, this is consistent with previ- 
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FIG. 7: I-V curves for BDT attached to Au contacts for different 
dss: (a) LDA and (b) ASIC results. In (c) and (d) we present LDA 
and ASIC calculated I-V curves where dss = 1.9 A but the angle 
between the molecule and the direction of the transport is 30°. 



ous results and it is due to the realignment of the HOMO 
of the molecule to of gold. Such a feature is investigated in 
more detail in figure[8]where we report the zero-bias T{E) for 
different dss- When dss is increased, the transmission peaks 
corresponding to molecular orbitals become narrower, as ex- 
pected, due to the weakened molecule-lead coupling. How- 
ever, charge transfer between the molecule and the metal is 
also affected, so that the charge on the actual molecule slightly 
increases as dss gets larger. This extra charge produces an 
electrostatic shift of the HOMO level towards the Au Fermi 
energy. As a consequence the narrowing of the resonance 
peak at the HOMO level is compensated by the upwards shift 
of the same peak, resulting in an enhancement of the zero-bias 
conductance. When ASIC is used this feature is strongly sup- 
pressed. In fact, the ASIC calculated transmission peaks are 
narrower and further away from that their LDA counter- 
parts. Thus, although also for ASIC the molecular HOMO re- 
aligns with respect to the Au Fermi level, this is is not enough 
to produce high transmission at reasonable dss- Note that a 
similar realignment is expected for weaker binding anchoring 
geometries which we will investigate in the remainder of this 
section (see figures l9b. [m[T2l and[T4ll 

The second contact geometry investigated is that where the 
S atom is connected to the bridge site on the fee (111) Au sur- 
face at both of the electrodes (figure [S}?). Total energy DFT 
calculations suggest that this configuration (dss = 2.09 A), 
has a lower energy than that of the hollow site^^ . We calculate 
a LDA zero bias conductance for the bridge site of 0.1 Go, 
which is lower than the value of 0.23 Go obtained for the hol- 
low site. This can be seen from the I-V curves presented in 
figure |9^. Interestingly this lower conductance found for the 
bridge site with respect to the hollow site is a low-bias feature 
and the two I-V curves matches closely for V > \.h Volt. The 
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FIG. 9: I-V curves for BDT attached to the Au (111) surfaces via 
a) bridge site and b) top site. For comparison in both cases the I-V 
curves for the hollow site geometry are also reported. 



FIG. 8: Transmission coefficients for BDT attached to Au contacts at 
different dss. The left plots correspond to LDA and the right ones to 
ASIC results. Note how the transmission peaks narrow and how the 
HOMO peak moves closer to the gold Ef as dss is increased. This 
realignment of the HOMO of the molecule has the effect of compen- 
sating for the weakening of the coupling, thus producing the counter- 
intuitive result of the low bias conductance increasing by enlarging 
dss- 



effect of ASIC on the transmission of the bridge site configu- 
ration is rather similar to that on the hollow site. Also in this 
case the molecular HOMO-LUMO gap opens and the current 
gets suppressed. The zero bias conductance for the bridge site 
is calculated to be 0.06 Go, the same value found for the hol- 
low site. Hence, whether the molecule is anchored to the hol- 
low site or to the bridge site makes relatively little difference 
to the low-bias transport properties as obtained with ASIC. 
This further confirms that our ASIC calculations are certainly 
more compatible than their LDA counterparts with the relative 
stability of the experimental conductance histograms. 

The top site of the Au fee (111) surface (figure [S]^) is the 
next anchoring site to be investigated. The equilibrium dis- 
tance dss is 2.39 A and the calculated I-V curves are pre- 
sented in figure |9l5. This time both LDA and ASIC present I- 
V characteristics with much higher conductance than that as- 
sociated with the hollow site. This is particularly dramatic for 
LDA for which the zero-bias conductance goes from 0.23 Gq 
for the hollow to 0.65 Go for the top site. For ASIC the in- 
crease is only twofold with the zero bias conductance going 
from 0.06 Go to 0.12 Go. Also for this case the general in- 
crease in current can be correlated with the larger Au-S bond- 
length of the top site compared to the hollow site. This has the 
same effect found when analyzing the hollow site under strain 
and it is due to the realignment of the HOMO of BDT to Ep 
of gold. It is important to point out that the top site is less 
energetically favorable than both the bridge and the hollow 
sites. Although it can still play a role in breaking junctions, 
we believe that this may be only of secondary importance and 
therefore would not have a considerable influence on the typ- 
ical conduction histograms. 

In actual breaking junctions, the electrodes are under strain 



and so the thiol groups are less likely to attach to a regular 
surface site. For instance in figure [TO] we present the case 
where one electrode connects to the molecule via the hollow 
site (dgg = 1.9 A) and the other via a gold adatom, with a 
distance of 2.39 A between the sulfur and the adatom. Re- 
cent X-ray standing wave expeiiments'"' demonstrate that gold 
adatoms may be the most favorable sites for the thiol groups in 
self-assembled molecular monolayer. This therefore may be a 
likely configuration for the bottom electrode in a STM-break 
junction. 




FIG. 10: (Color on line) BDT molecule attached asymmetrically to 
two Au (111) surfaces. On one side the bonding is at a hollow posi- 
tion (dss = 1.9 A), while on the other is to a gold adatom (sulfur- 
adatom separation=2.39 A). Color code: Au=yellow (light gray), 
C=black, S=dark yellow (gray), H=blue (dark gray). 

The orbital resolved DOS, transmission coefficients and I- 
V curves for this system are presented in figure [TT] This con- 
figuration shows the largest difference between the conduc- 
tance calculated with LDA and that calculated with ASIC. The 
LDA conductance is 0.32 Gq, whereas when ASIC is applied 
this drops by a order of magnitude to 0.03 Gq. Because of the 
weaker interaction between the tt orbitals of the molecule and 
the gold surface, the molecular orbitals in the DOS (figurefTTh 
and figure [TTb) and hence the peaks in the transmission coef- 
ficients (figure [TTt and figure [TTtl) are considerably narrower 
than for the case of the hollow site. These narrow levels are 
however closer to E-p, resulting in the relatively high conduc- 
tance at low bias. This originates from the same mechanism 
producing a change in conductance when dss is increased (see 
figure[8]). At variance with the case of stretched hollow bond- 
ing geometry, this time the ASIC does not shift the conduc- 
tance peaks enough to bring them close to the Fermi level of 
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gold, and the conductance remains rather small. 
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FIG. 11: Transport properties of a BDT molecule attached asymmet- 
rically to two Au (III) surfaces according to the anchoring geometry 
of figure [To] The left plots correspond to LDA and the right ones to 
ASIC. The upper panels are the DOS of the S and C tt orbitals ((a) 
and (b)), the middle are the transmission coefficients as a function 
of energy for various bias ((c) and (d)) and the lower are the I-V 
curves. Figure (f) is a zoom of (e) and compares our results with 
experiments from reference** . The vertical lines in (c) and (d) mark 
the bias window. 

As the peaks in the transmission are rather naiTow we 
checked our results against the choice of basis set. In figure 
[12] we present the orbital resolved DOS, transmission coeffi- 
cients and I-V curves for the enlarged 5d6s6p basis, to be put 
in perspective with those calculated with 6 s only in figure [TT] 
The most notable difference is the appearance of transmission 
away from E-p, in particular at low energies. This is connected 
to the gold d electrons, which are neglected in the small ba- 
sis. Importantly however the low-energy features around E-p 
are only slightly affected with the zero-bias conductance as- 
suming values of 0.47 Gq and 0.05 Go for LDA and ASIC 
respectively. This confirms the feasibility of using a restricted 
basis set for the I-V curve at low bias. 

Finally, we investigate two more configurations. The first 
considers adatoms as the bonding site at both sides of the 
junction (figure[T3]l, while the second investigates hollow sites 
where the two hydrogens of the thiol groups are not dissoci- 
ated in forming thiolate as in all the other cases studied (figure 
[TSl l. In the first case we expect a rather weak bond, which is 
confirmed by the sharp peaks in transmission shown in fig- 
ure [14] The LDA low-bias conductance, calculated with the 
5d6s6p basis for gold, is 0.43 Go, which drops to 0.19 Go 
when ASIC is applied. In this case the molecular HOMO is 
pinned at Ep and the transmission remains large. Interestingly 
this is a situation similar to the one we discussed in our pre- 
vious work in the context of tight-binding Hamiltonian-^. In 
this case the suppression of the conductance may further orig- 
inate from the derivative discontinuity of the potential, which 








30.0 



5-4-3-2-101234 5-5 -4-3-2-1012345 
(eV) (eV) 

20.0 r 






f) 


1 








LDAspd 

-- ASICspd 
Exp. 

















0.2 0.4 0.6 

V (Volts) 



FIG. 12: Transport properties of a BDT molecule attached asymmet- 
rically to two Au (1 1 1) surfaces according to the anchoring geometry 
of figure[TO] Results are presented for calculations using an enlarged 
5d6s6p basis for Au. The left plots correspond to LDA and the right 
ones to ASIC. The upper panels are the DOS of the S and C tt or- 
bitals ((a) and (b)), the middle are the transmission coefficients as a 
function of energy for various bias ((c) and (d)) and the lower are the 
I-V curves. Figure (f) is a zoom of (e) and compares our results with 
experiments from reference^. The vertical lines in (c) and (d) mark 
the bias window. 



may move the peak in the conductance away from Ep. Since 
this is not described by the ASIC scheme, we believe that the 
cuiTent in this case may be erroneously large. 





FIG. 13: (Color on line) BDT molecule attached to adatoms at both 
of the Au (III) surfaces. Color code: Au=yellow (or light gray), 
C=black, S=dark yellow (or gray), H=blue (or dark gray). 

The results for the case where BDT is attached to the hol- 
low sites but hydrogen atoms remain attached to thiol anchor- 
ing groups are presented in figure [16] The total energy for 
this situation is 1.47 eV higher than that obtained when the 
end groups are thiolate and the remaining two H atoms form a 
H2 molecule. As can be seen from the DOS in panels (a) and 
(b) and the transmission coefficients in panels (c) and (d), the 
transport is now through the LUMO of the system, the energy 
of which is lowered only slightly by ASIC. Hence, the ASIC 
zero-bias conductance of 0.09 Gq is higher than the value of 
0.06 Go calculated with LDA. However it is important to note 
that the small, but significant, downshift of the LUMO state 
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FIG. 14: Transmission coefficients for a BDT molecule attached to 
adatoms on both the gold (111) surfaces, calculated with the 5d6s6p 
basis for gold. Note that the HOMO is now pinned at the gold Ef . 




FIG. 15: (Color on line) BDT molecule with H atoms attached to the 
S atoms in the thiol groups (dss = 1.9 A). Color code: Au=yellow 
(or light gray), C=black, S=dark yellow (or gray), H=blue (or dark 
gray). 



caused by ASIC is only an artifact of our atomic approxima- 
tion to SIC. As an empty state, the SIC for the LUMO is not 
defined and no corrections should be applied. 



B. Benzenedimethanethiol 

The second device studied consists of a ben- 
zenedimethanethiol (BDMT) molecule attached to two 
gold fee (111) surfaces. Two different isomers of are inves- 
tigated; in the first the sulfur atoms of the thiol groups lie in 
the same plane of the benzene ring (figure [TTk). while in the 
second they are aligned out of plane (figure [T7b). In both 
cases they anchor the molecule to the hollow site of the gold 
fee (111) surface. 

The orbital resolved DOS, transmission coefficients and I- 
V curves for the first isomer are presented in figure[T8]for both 
LDA and ASIC. In this case, as one can see from the orbital 
resolved density of states (figures [TSb andfTSb). the HOMO- 
LUMO gap is much larger than that of BDT. The ASIC again 
has the effect of shifting the occupied orbitals downwards in 
energy, but this has little effect on the transport since the res- 
onances in the transmission coefficients due the HOMO and 



O 
Q 



1 1 II II 1 1 1 1 1 

a)jl 




— S 71 

-- C Jt 


1 J 








■ 










^ ■ 




™ ■ 





M 



-5 -4 -3 



^40 
^ 0, 



■2-10 12 
E-E^ (eV)r 



3 4 5-5 -4 -3 -2 -1 1 2 3 4 5 



LDAH 

— - ASICH 
-.. LDA 
. . - . ASIC 



< 







0.5 1 1.5 

V (Volts) 



(eV) 


LDAH 

— - ASICH 
-.. LDA 
.-. ASIC 
Exp. 


:f) 





0.2 0.4 0.6 

V (Volts) 



FIG. 16: Transport properties of a BDT molecule attached to two 
hollow sites at both sides of the junctions. In this case the H atoms 
of the thiol groups are left bounded to S and do not dissociate in 
forming a thiolate end group. The left plots correspond to LDA and 
the right ones to ASIC. The upper panels are the DOS of the S and C 
TT orbitals ((a) and (b)), the middle are the transmission coefficients 
as a function of energy for various bias ((c) and (d)) and the lower 
are the I-V curves for the system both with and without the hydrogen 
atoms attached to the sulfur. Figure (f) is a zoom of (e) and compares 
our results with experiments from reference^. The vertical lines in (c) 
and (d) mark the bias window. 



FIG. 17: BDMT isomers attached to the hollow site of the Au (1 1 1) 
surface. The sulfur- surface distance is 1.9A. Color code: Au=yellow 
(or light gray), C=black, S=dark yellow (or gray), H=blue (or dark 
gray). 



LUMO lie well outside the bias region investigated. 

Therefore, ASIC increases the HOMO-LUMO gap but the 
actual transmission coefficient around E-p does not change 
much from its LDA value. The ASIC conductance at zero bias 
is calculated to be 0.004G'o, compared to O.OO6G0 obtained 
with LDA only. Both these values are one order of magnitude 
larger than the experimental value of O.OOO6G0 obtained by 
Xiao et. al.^. It also follows that the I-V curves (figures [TSfe 
andfTSlF) obtained with LDA and ASIC respectively are quite 
similar. 

The I-V curve calculated for the second isomer is presented 
in figure[T9l Also in this case the effects of ASIC are minimal 
due to the large HOMO-LUMO gap. The conductance at zero 
bias is calculated to be 0.0 15 Go with LDA, compared to the 
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C. Biphenyldithiol 

The third and final system investigated is that of 
biphenyldithiol (BPD) attached to two gold surfaces as shown 
in figure l20l Once again we consider only the gold fee (111) 
hollow site as the anchoring site. The planes of the two ben- 
zene rings are usually tilted with respect to each other by an 
angle known as the torsion angle, that in this case is calculated 
to be 37°. 




FIG. 18: Transport properties of a BDMT molecule attached to the 
gold (111) hollow site for the isomer of figure [T7] The left plots cor- 
respond to LDA results and the right ones to ASIC. The upper panels 
are the DOS of the S and C tt orbitals ((a) and (b)), the middle are the 
transmission coefficients as a function of energy for various bias ((c) 
and (d)) and the lower are the I-V curves. Figure (f) is a zoom of 
(e) and compares our results with experiments from reference**. The 
vertical lines in (c) and (d) mark the bias window. 



value of O.OlSGo for ASIC. Interestingly this is about a factor 
3 larger than that calculated for the first isomer. 
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FIG. 19: I-V curves for both BDMT isomers attached to gold con- 
tacts. In both cases ASIC has only marginal effects on the currents, 
which differ by a factor of about three for the two different isomers. 



In general, in contrast to the case of BDT, for BDMT ASIC 
does not appear to improve the agreement between theory and 
experiments. The fact that the Fermi level of Au is pinned in 
the middle of the large HOMO-LUMO gap of BDMT, which 
is already well described by LDA, makes the ASIC correc- 
tions rather marginal to transport at moderate voltages. Given 
the fact that the geometry of the anchoring is not known and 
that we have only investigated one case, we may conclude that 
the disagreement between theory and experiments may simply 
lie in the details of the anchoring geometry. 



FIG. 20: BPD molecule attached to the hollow site of the Au (HI) 
surface. The sulfur- surface distance is I.9A. Color code: Au=yellow 
(or light gray), C=black, S=dark yellow (or gray), H=blue (or dark 
gray). 

The orbital resolved DOS, transmission coefficients and /- 
V curves for this system are presented in figure |2T| for both 
LDA and ASIC. Once again, ASIC has the effect of lowering 
the energy of the occupied molecular orbitals, as can be seen 
from the DOS (panels (a) and (b)). This in turn results in open- 
ing up the conductance gap in the transmission as shown in 
panels (c) and (d) of figure|2Tl For this molecule, the HOMO 
is near as in the case of BDT, giving a LDA conductance 
of 0.07Go at zero bias. When ASIC is appHed, the HOMO 
is shifted downwards in energy and outside the bias window 
for voltages up to 2 Volt. In this case the ASIC conductance 
drops to O.OI8G0 at and the whole low bias I-V curve 
shows rather smaller current with respect to its LDA counter- 
part. 

The optimum torsion angle is 37°^^. However, this may 
fluctuate due to temperature or when the molecule is under 
strain in a breaking junction. Panel (e) of figure |2T| shows 
the 1-V curves calculated for the equilibrium torsion angle of 
37°, whereas panel if) shows the result for the case when the 
benzene rings are in the same plane (i.e. when the torsion 
angle is 0°). Reducing the torsion angle causes an increase 
in the transmission since the overlap between the tt orbitals 
is increased. The conductance at zero bias for a torsion angle 
of 0° is 0.09Go with LDA only, and 0.024Go when ASIC is 
applied. 

These results show that ASIC has an effect on BPD similar 
to the one it has on BDT, i.e. it shifts the HOMO downwards 
in energy and reduces the zero-bias conductance. In this case 
however even our ASIC results differ from the experimental 
data of Dadosh et. ali^S by several orders of magnitude. Since 
we are essentially in a tunneling situation, the magnitude of 
the current is severely dependent on the tunneling matrix ele- 
ments, which in turn are rather sensitive to the details of the 
anchoring geometry. Although the hollow site is believed to 
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FIG. 21: Transport properties of a BPD molecule attached to the 
gold (111) hollow site. The left plots correspond to LDA and the 
right ones to ASIC. The upper panels are the DOS of the S and C vr 
orbitals ((a) and (b)), the middle are the transmission coefficients as 
a function of energy for various bias ((c) and (d)) and the lower are 
the I-V curves. Figure (e) is the I-V curve for a torsion angle of 37° 
and (f) is for a torsion angle of 0° . The vertical lines in (c) and (d) 
mark the bias window. 



be the preferential site, such a large discrepancy between the- 
ory and experiments could be attributed to the contact geom- 
etry. 



IV. CONCLUSION 

In view of the many contradictory results for the transport 
properties of molecular junctions, first principles computa- 
tional tools are becoming increasingly important. These, how- 
ever, also suffer from some uncertainty in the results, rooted in 
the approximations used for computing the junction electronic 
structure. In this work, we have addressed the problem of the 
self-interaction error in DFT and we demonstrate how a sim- 
ple self-interaction correction scheme can be used to improve 
the agreement between theory and experiments. In particu- 
lar, we have considered the problem of correctly reproducing 
the molecular ionization potential as the negative of the single 
particle Kohn-Sham HOMO energy. 

Our ASIC scheme is capable of achieving such a descrip- 
tion, and for most of the molecules investigated this has the ef- 
fect of moving conductance resonances away from the Fermi 
level of the electrodes. Thus, molecular junctions with zero 
bias resonant transport through the HOMO become insulating 
when calculated with ASIC. We then investigated transport 
through BDT attached to the gold (111) surfaces in several 
different anchoring geometries. In general, we have found a 



systematic reduction of the low bias current when ASIC is ap- 
plied. Most importantly, we have demonstrated that a number 
of anchoring geometries all yield similar low bias conductivi- 
ties, which in most cases are rather insensitive to bond stretch- 
ing or bond angle changes. This suggests that several anchor- 
ing geometries might contribute to the peaks in the typical ex- 
perimental conductance histograms, thus providing a natural 
explanation for their stability. 

Nevertheless, even with ASIC, some problems still remain. 
In general, if the first molecular orbital to enter the bias win- 
dow as the bias increases is the LUMO, then our method po- 
tentially fails. In fact, ASIC only partially restores the DFT 
deiwative discontinuity, and erroneously Ehomo 

(N). Since the size of the derivative discontinuity is 



"LUMO 



unknown, it is difficult at this time to evaluate the impact of 
such an error. A second problem is connected to the accuracy 
of the calculation of the wave-functions. When no molecu- 
lar levels appear at the Fermi level, the transport is essentially 
tunneling like, and therefore becomes crucially dependent on 
the accuracy of the calculated hopping integrals between the 
molecule and the electrodes. ASIC does not drastically im- 
prove the quality of the LDA wave-function and thus one 
might expect some eiTors in the tunneling transmission coeffi- 
cients. Additionally, we do not correct the surface atoms of the 
electrodes for SI, so it is expected that our ASIC results in gen- 
eral produce overextended orbitals and hence a systematically 
overestimated current. Finally, it is worth noting that ASIC 
still overestimates the polarizability of molecule s.22il£, with a 
quantitatively incorrect prediction of the response exchange 
and correlation field. The use of XC potentials constructed 
from exact charge densities'*' which coiTect both the molecule 
and the metallic surfaces may offer a solution to these prob- 
lems. However, for some of the molecules investigated the 
disagreement between theory and experiments is still of sev- 
eral orders of magnitude. This is difficult to associate with 
any systematic eiTor in the theory, and might be connected to 
the specific contact geometry. 

In conclusion, NEGF combined with the ASIC electronic 
structure method offers a description of transport through 
molecules which is considerably improved with respect to 
standard NEGF LDA/GGA. Since this is also a computation- 
ally undemanding scheme, it can be employed for large sys- 
tems, and therefore can become an important tool for system- 
atically predicting the performance of molecular junctions. 
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